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^ ■ Abstract 

o; 

' 55 ■ We analyze via theoretical approaches and molecular dynamics simulations the collective mode 

structure of strongly coupled two-dimensional binary Yukawa systems, for selected density, mass 
and charge ratios, both in the liquid and crystalline solid phases. Theoretically, the liquid phase is 
described through the Quasi-Localized Charge Approximation (QLCA) approach, while in the crys- 
. talline phase we study the centered honeycomb and the staggered rectangular crystal structures 

\Q ' through the standard harmonic phonon approximation. We identify "longitudinal" and "trans- 



verse" acoustic and optic modes and find that the longitudinal acoustic mode evolves from its 
weakly coupled counterpart in a discontinuous non-perturbative fashion. The low frequency acous- 
tic excitations are governed by the oscillation frequency of the average atom, while the high fre- 



quency optic excitation frequencies are related to the Einstein frequencies of the systems. 
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I. INTRODUCTION 



Yukawa systems, i.e. many particle systems where the pair interaction potential energy 



is 



<P{r) 
(p(r) 



ZxZ 2 ip{r) 
exp(— Kr) 



(1) 



have been of interest for some time. The Yukawa potential has the unique feature that 
by varying the screening parameter k the potential can assume the feature both of a short 
range (hard sphere like) and of a long range (Coulomb like) interaction potential. Since 
the 1970-s this feature motivated a number of investigations relating to the properties of 
Yukawa liquids and solids, their phases and phase transitions [l ^. Quite apart from this 
academic interest, the Yukawa potential has been recognized as a good approximation for the 
interaction potential between charged particles in colloids [4] and, more recently, in complex 
(dusty) plasmas, where the original Coulomb interaction between the main constituents 
is transformed by Debye screening into a Yukawa-type potential (for a review see, e.g. 
. Complex plasmas constitute an especially suitable medium for the study of waves 



and collective excitations because these are much less damped than in colloidal systems. In 



recent years a host of papers, both theoretical 



collective modes in Yukawa systems 



15] and experimental have studied 



18H25I]. Most of the experimental work on colloidal 
systems and complex plasmas has focused attention on two dimensional (2D) layers. We 
note that the physics of the 2D and 3D systems, the dynamics of the collective excitations 
in particular, is, in fact, quite different and addressing 2D and 3D systems separately is also 
warranted on theoretical grounds [4|. 

The strength of the coupling between the particles can be characterized by the nominal 
bare coupling parameter T = (Z 2 e 2 )/(k^Ta), with a being the Wigner-Seitz radius. A 
physically more meaningful T e g that, basically represents the ratio of the potential and 
kinetic energies, can be defined for orientation purposes as r e fj = T exp(—na) 2^|, although 
more sophisticated expressions are available 27H29|. 

The main interest lies in the behavior of the strongly coupled state, r e fr ^> 1. In this 
strongly coupled state the system can be either in the dense liquid or in the crystalline solid 
phase. 



Past works on the dynamics have overwhelmingly concentrated on Yukawa systems con- 
sisting of one single component, the equivalent of the One Component Plasma, both in 
three and two dimensions (Y3dOCP and Y2dOCP, respectively). A great deal of theoreti- 



nn 



10 



32 



- |35| efforts have been devoted to the 



cal [9|, lllMl5l. l30l . l3lj and computer simulation 
mapping and understanding of the collective mode structures in these systems, both in the 
liquid and solid phases. The theoretical methods required in the two situations are, of course, 
quite different. Once the lattice structure is identified, the crystalline solid is amenable to 
the standard harmonic phonon analysis. Concerning the treatment of the collective modes 
in the strongly coupled liquid phase, the Quasi Localized Charge Approximation (QLCA) 



approach developed by Kalman and Golden [36j, |37j has turned out quite successful. The 



QLCA borrows ideas from the harmonic phonon approximation and describes the liquid in 
terms of particles trapped and oscillating in local potential wells of the fluctuating poten- 
tial (for details, see |37|). As a result of these works, the collective mode spectra of the 
Y3dOCP and Y2dOCP are well understood and this understanding is well corroborated by 



observations 



16 



As to strongly coupled Yukawa mixtures consisting of more than one single species, in 



ncsc 



particular binary Yukawa mixtures (Y3dBM, Y2dBM), the collective dynamics of t 
systems constitutes a largely unexplored area (see, however, a recent work by Daligault 
, even thou gh t he problem is of great theoretical interest. (For a related one dimensional 



3* 



problem see 



39 



40]). One expects that the simple analytic structure of the Yukawa potential 



allows one to derive nearly exact solutions, which will elucidate the common features of the 
dynamics of binary liquids and solids 



41 



42j. Also, the flexibility of the Yukawa interaction 



would make the qualitative features of the results to serve as paradigms for the collective 
mode structures of binary systems interacting through other potentials as well (alloys, dipole 
systems, etc.). From the point of view of actual applications, the creations of binary complex 
plasmas has technical problems, but, nevertheless, one expects that such strongly coupled 
complex plasmas of two different grain species will become available in the near future. 

This paper, the first in a series, presents a systematic study of he collective mode spectra 
of the Y2dBM system. The system consists of two species, with charges Z\ and Z2, masses 
mi and 1712 and densities rt\ and 712 (or concentrations c\ and C2), respectively. Our strategy 
is similar to the one followed in our previous works on the Y3dOCP and Y2dOCP: for the 
theoretical analysis of the liquid state we apply the QLCA formalism; for the crystalline 



solid we calculate dispersion relations by the standard method. In both cases, we parallel 
our theoretical analysis with detailed Molecular Dynamics (MD) simulations of the den- 
sity and current fluctuation spectra of the system; it is, then, the positions of the peaks 
of the fluctuation spectra from which the dispersion relations are inferred. It has to be 
noted, however, that following this road map is fraught with questions stemming from the 
fundamental difference between the binary and single component systems. Concerning the 
QLCA, is it justified to represent the system through separate collective coordinates for 
each of the species in a liquid where the two species are spatially not separated? Concerning 
the MD, if different partial fluctuation spectra provide conflicting information, which one of 
them should be accepted as most relevant to the actual dispersion? Finally, one has to be 
aware of the fact that in the presence of different Z\ and Z 2 charges with different C\ and c 2 
concentrations, the liquid phase is governed by a complex phase diagram 4ll447] in which 



only certain combinations of these parameters allows a homogeneous system. In the solid 
)hase, similarly with a given set of parameters only certain lattice structures are permissible 



48 



49|. 



We tackle these issues along the work as presented in the sequel. We have to empha- 
size though, that our goal is restricted to determining the existence, interrelationships and 
dispersion of the collective modes. We do not address a number of related problems: the 
damping of the modes, the detailed structures and the link between the various fluctuation 
spectra, the critical freezing values of T, the nature of the underlying order in the liquid 
phase, lattice stability and structures, etc. The issues investigated in this paper are orga- 
nized according to the following plan: Section II is devoted to the description of the liquid 
phase and Section III of the crystalline solid phase. In each case we first analyze the qual- 
itative features of the optic (u finite at k = 0) and then the acoustic (u — > as k — > 0) 
excitations, before presenting the description of the full mode structure. In Section III we 
compare the mode structures in the two phases and draw conclusions. (For a preliminary 



account of some of the results pertaining to the optic modes see [50] and to the acoustic 



modes see |51|). 



Whenever not noted otherwise, we measure frequencies in units of u\, the plasma fre- 
quency of species 1, use Ti, the bare coupling value for species 1 to characterize the coupling 
strength, and adopt na — 1 (a = ^aia 2 ) for the screening parameter. 



II. STRONGLY COUPLED LIQUID PHASE 



The theoretical analysis of the mode structure in the liquid state is based on the QLCA 
approach, as discussed above. The fundamental equation for the dynamical matrix is 

2 \fn A n B 



C^(k) = -Z A Z B e< 



y/m A m B 
Z c n c 



j d 2 r|^(r) (exp(-ik • r) - 5 AB ) [1 + h AB {r)\ - 



6 ** E y^" v{j) [1 + hAc{r)] 



L J d 2 f^(r) (exp(-ik • r) - 5 AB ) [1 + h AB {f)\ + 
6 AB &AC j d 2 f ^(r) [1 + h AC {f)\ (2) 



with 



W»{r) = d^dMr), (3) 
where ip(r) is the Yukawa interaction, <p(r) = exp(— nr)/r characterized by the screening 



constant k. Then 



a{y) = l + y + ^y 2 , b(y) = l + y. (4) 

Additional notational conventions are 

2 _ 2ne 2 Z A Z B n B 
m A a 

2 2ixe 2 Z A Z B ^/n A n B 

u AB — , 

^/m A m B a 

a = y/aia 2 
a A = l/y/nn A 
R = Kd 
f = r/a 

y — Kr (5) 

with Z, m, n and a representing the charge number, mass, density and Wigner-Seitz radius 
for the respective components. The Q AB and u AB frequencies are the nominal Einstein and 
nominal plasma frequencies of the system. The h AB pair correlation functions are to be 
obtained from the MD simulations, as described below. 



The elements of the C-matrix can be expressed in terms of the kernel functions /C, C: 
C L AB = co 2 AB J ^)C(kr, y) [1 + h AB (f)\ - 

5 AB ^Icj %^y) [l + hacir)] 



C(all) 

rp£(kr,y) [1 + h AB (f)\ - 
^ab E tf BC / J^(0, y) [1 + h BC {f)\ (6) 

C(all) J T 

with the kernel functions given by 

}C(u, r) = - exp(-y) { [l + y + y 2 } J (u) - 3 [l + y + y 2 /3] J 2 (u) } 

C(u, r) = - exp(-y) { [l + y + y 2 ] J (u) + 3 [l + y + y 2 /3] J 2 (u) } (7) 

In order to clearly display the behavior in the vicinity of k = we also introduce 

Q(u, r) = K.(u, r) — /C(0, r) 
H(u,r) = C(u,r) - C(0,r) 

F(r) = -/C(0,r) = -£(0,r) (8) 

The integrals of the kernel functions over the pair correlation functions 1 + h(r) are 

K AB (k) = J ^/C(£x,r) [1 + h AB (r)\ 

L AB (k) = J ^C(kr,r) [1 + h AB (r)\ 

F ab = J ^Hr)[l + h AB (r)} (9) 

These integrals would be divergent at r = 0, were it not for the pair correlation function 
1 + h(r) that becomes at r = 0. Similarly 

G AB (k) = J ^G(kr,r) [l + h AB (r)\ 

H AB (k) = J ^n(kr, r) [1 + h AB (r)\ (10) 
Introducing the asymmetry parameters p and q 

p 2 = Z 2 n 2 /Z 1 n 1 

q 2 = Z 2 m 1 /Z 1 m 2 (11) 



one obtains for the longitudinal elements 

^nW = y [G n (A:)+p 2 F 12 ] 

C 2 L 2 (k) = ^-[p 2 q 2 G 22 (k) + q 2 F 12 ] (12) 
while the transverse elements are 

CUk) = ^[H ll (k) +P ^F 12 ] 
C T l2 {k) = ^-pq [H 12 (k)-F 12 ] 

CUk) = ^[pVH 22 (k) + q 2 F 12 ] (13) 

We have found it useful to introduce u>i (= Wn) as the reference frequency. In general, 
there exist 4 modes as the roots of the characteristic equations, 

\\C L A £-co 2 \\=0, (14) 

which will be labeled u;+, u;+, uit, oo 1 -- The ± notation identifies the polarizations in species 
space of the modes: the "+" sign designates polarization where the two components move 
in-phase, while the "-" sign designates polarization where the two components move out- 
of-phase. The two + modes are acoustic (uo — > as k — > 0) and the two - modes are 
optic modes (u finite for k — 0). In addition, the modes are labeled as Longitudinal L or 
Transverse T, referring to the their polarization with respect to k when the propagation is 
along the principal axes. We note that the elements of the C-matrix, and consequently the 
eigenfrequencies, depend only on the two p and q combinations of the originally introduced 
three Z — Z 2 /Z 1 , M — m 2 /m 1 , N — n 2 jn\ parameters. 



A. Optic modes 

At k = 0, GAB(k) oc HAB(k) oc 0(k 2 ), thus u~(k = 0), the gap frequency, is longitudi- 
nal/transverse degenerate, as it should be for an isotropic liquid: 



co GAP = ^ = 0)= ui(k = 0) = coi\/-(p 2 + q 2 ) F] 



1 



'12 



= \j\ (n 2 l2 + n 2 2l ) f 12 = {n 2 2 + ti 2 21 ). (is) 



In view of Eqs. §6§ through (jUJ) Fab can be interpreted as the average potential generated 
by species B in the environment of a particle of species A. 



Fab = 1 LJ d2f <*(r))[l + Wr)] 



O 2 



(16) 



with (...) designating angular averaging. The Qab frequency represents the oscillation 
frequency of a particle of species A in the frozen environment of particles of species B. We 
note that it is the correlation dependent Q-s, rather than the nominal Q-s that are the real 
Einstein frequencies of the system 34(, with a similar definition being used in the theory 



of liquids 



521 ] . In a single component system the Einstein frequency Q also provides the 



cu(k — > oo) limiting frequency 11]. 

In order to find the k — > oo limits for the binary systems we re-express the elements of 
the C-matrix as 



°22 



IjJ 



12 



{K ll {k) + {n\ l + ^\ 2 )} 
K 12 (k) 

{K 22 (k) + (n 2 22 + n 2 21 )} 



(17) 



In the k — > oo limits the .fT-terms vanish. This can be seen by observing that 

/j 
-^)C(u,y)[l + h AB (u/k)} (18) 

and that [1 + h(r — > 0)] — > fast enough to make this happen. Similar considerations apply 
to the transverse elements. 

Thus the k — > oo respective upper and lower effective Einstein frequencies become fli, flu: 



uj-{k — > oo) 
uj + {k — > oo) 

The calculated gap frequencies and the effective Einstein frequencies as functions of T, 
together with the results obtained by MD simulations (see below) are shown in FIG [TJ also 
shown is the variation of the correlation integral F\ 2 . In anticipation of the results of the 
next Section, we have also indicated the gap frequencies in the crystal lattices. We will 
further comment on the relationships between these gap frequencies in the next Section. 

8 



~ (nf ! + ny = Cij 



(19) 
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FIG. 1. Liquid state: QLCA gap (•) and Einstein (■) frequencies vs. F. The arrows indicate 
the positions of the corresponding gap in the lattice. The Inset shows the V dependence of the 
correlation integral F12 (symbols), which does not vary with the mass ratio, (a) ri2 = n\/2; (b) 
n 2 = n\. 



B. Acoustic modes and sound speed 



We now turn to the calculation of the acoustic modes in the binary system. We are 
interested primarily in the small-/c behavior, which will lead to the determination of the 
sound speed. 

First we observe that by dropping h(r) in the integrals for the G(k) and H(k) functions 



9 



the resulting G°(k) and H°(k) integrals become doable and provide the RPA expressions: 



u 



- exp(-Kr) [{l + y + y 2 } {1 - J (u)} - 3 {l + y + y 2 /3} J 2 („ (J 

r VK + fir 

^°( fc ) = / ^ex V (~Kr)[{l + y + y 2 }{l-J (u)} + 3{l + y + y 2 /3}J 2 (u)} =0 (20) 



The C-matrix equivalent to the cold RPA approximation would be obtained by dropping 
the F12 terms in (TT2"]) and (1T3"]) . and using ( 12"U|) for Cn, Ci 2 and C 2 2- Then one obtains the 
RPA result 



L k 
UJ, = u - 



Vk 2 + k 2 



LU L =0 



oo = wi VI + P 2 ? 2 = V^i 2 + w 2 2 - (21) 



Note that the intuitively more reasonable requirement that in order to obtain the RPA 
limit one sets hn equal to zero everywhere in ( TT2l) and ( FTBl would result in a meaningless 
divergent integral for F\2- This feature shows that there is no smooth transition from the 
QLCA expression to the RPA. In other words, in contrast to the case of the YOCP, in 
the YBM the RPA Eq. fj20|) cannot be simply amended by adding correlational corrections 
in order to obtain the strong coupling expression: the strong correlations show up in an 
essentially non-perturbative fashion. 

Returning now to Eqs. (fl2"j) and f fl3|) we now calculate the small k expansion. The result 
is given in terms of the integrals 



5 f°° 

Uab = 'V,L dy 



3 2 

l + y + -y 2 



exp(~y)h AB (r), 

v AB = —^J o d v[ 1 + y- y 2 } ew(-y)h AB (r). (22) 

Thus the longitudinal and transverse C\ B and C T AB matrix elements in the k — > limit 



10 



become 



C&(Jfc -> 
<7&(* 
C&(A; -> 
C£(A: ->• 
C&(A; -> 
C T 22 {k -> 



^{ (1 _ f/n) | + ^ Fl2 } 

w 2 f A; 2 1 1 

y S P9(l - ^12) -3- - ^pq^F 12 > 

f{ P V(i-fe)^M 2 } 



y n — +p z rf 12 

2 K 



~2 



2 r k 2 1 \ 
}\pqV 12 — - -pqKF 12 \ 



(23) 



Proceeding now from (1231) . after some algebra one finds the small- k expansion of the 
relevant u>+(k), u+(k) mode frequencies as 



(^) 2 (A;->0) =w 2 |l 



U n + 2p 2 U 12 + p 4 U 22 ] k 



(l+p 2 Y 



H 



(ioi) 2 (k -» 0) = 



_ 2 (Vn + 2p 2 V 12 +p 4 V 22 \ k 



(24) 



i+p 2 r j « 

While the first term in (cuf ) 2 is RPA-like in appearance since it shows no explicit depen- 
dence on h(r), in fact it reflects an essentially strong coupling behavior, the correlational 
effects manifesting themselves through the Q coefficient, which we will refer to as the "virtual 
average atom" (VAA) frequency (this f req uency has also been mentioned in relation to the 
self-diffusion coefficient of a plasma in [53] ) . 

,2 



-2 2 1 

U = LO 



2\ 2 



1 p 2 + q 2 



(1+p 2 ) 



(25) 



The Virtual Atom in fact represents an entity created from the averages of the system 
parameters. To see this, Eq. [23 is re-written in terms of the average charge and mass as 



00 



! 2ne 2 (Z) 2 



a (m) 
n = rii + n 2 . 



■n. 



The averages are defined through 



(X) 



Ysj X i n i 



(26) 



(27) 
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Compare now u with of Eq. [2TJ the dramatic difference in the dependence on the plasma 
parameters, in particular on the mass ratio, is evident. (A similar result but restricted to 
the Zi = Z 2 case was already anticipated in 51]). 

The notion of the VAA originates from the literature, pertaining to liquid alloys and dis- 



ordered binary systems 



54H56| . as a heuristic concept. Here the derivation of this behavior, 



as a result of the evolution of the system from weak to strong coupling, is given. 

All the observations now made on the k — > behavior of the acoustic mode can be 
translated into statements about the sound speeds 



uj+ T {k) /k 



J fc->0 



(28) 



Thus, according to (|2"3|) and (12"5|) . the longitudinal sound speed at weak coupling has its 
RPA value, governed by Uq] for strong correlations the sound speed is substantially reduced 
and strong correlations manifest themselves, in contrast to the YOCP, in two ways: first, 
by morphing the mean field contribution into one whose properties are dictated by the VAA 
and do not explicitly depend on the correlations and, second, by generating an explicit 
correlational correction. For the transverse sound speed, similarly to the YOCP, there is no 
/i-independent contribution. 

In parentheses we remark that to what extent the weak coupling value of the sound 
velocity is well represented by the RPA (or "cold fluid" ) expression is not clear. It is generally 
assumed that it is j^, 0|. Nevertheless, the issue is that while for a Coulomb system there 
exists a clear rigorous derivation (also supported by ample observational evidence) that shows 
that in the r — > limit the RPA is correct, no such demonstration is currently available for 

58j that for a finite range system the 



a Yukawa system. In fact, there is reason to believe 
description of the behavior of the system in the weak coupling limit is more involved. All 
this, however, has very little bearing on our conclusion that the sound speeds and the low 
frequency excitations in the strongly coupled system are governed by the frequency of the 
VAA and thus are quite different from their weak coupling counterparts. 

We have studied the T-dependence of the sound speeds and of the related effective masses, 
the latter being defined by subtracting the explicitly correlation dependent term from the 
sound speed coefficient 

m cS ufa 2 (Z) 2 



mi 



n 2 

1 + — 

m 



(l-U) 



(29) 
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by MD simulations for the parameter set given previously. Results are shown in Figs. [2] and 
[3] for T values between T = 5 and T = 120. For T=l and T=5 sound speed values calculated 
through the (Vlasov Equation based) RPA approach are also displayed. 
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FIG. 2. Liquid state: Longitudinal sound speeds. (•) MD, line: QLCA, gray shaded area/line: 
lattice value. For T = 1 and T = 5 the RPA (Vlasov equation based) values of the sound speeds 
are also indicated (■). (a) ri2 = ni/2; (b) n-2 = n±. 



At the high T end the QLCA predicted behavior is in excellent agreement with simulation 
results. As T approaches the freezing boundary, the sound speeds smoothly join their values 
in the crystal lattice, which are also given, in anticipation of the results of the next Section. 
Some further comments on the relationship between the sound speeds in the two domains 
will be given there. In the liquid, as T is lowered, the remarkable decrease of the effective 
mass and the concomitant increase of the sound speed can be observed. At the same time, 
the QLCA sound speed, in general, stays below the observed value because the QLCA 
ignores the modification of the effective mass as T is reduced. It can be noted, that even at 

13 



the relatively low T = 5 value the strong coupling behavior seems to be still dominant and 
the sound speed is much below its calculated RPA value. The behavior of the sound speed 
below this T value is not clear: it is a domain that would require substantial theoretical, 
simulation and experimental work to arrive at a reliable and coherent picture. 




FIG. 3. Liquid state: Effective masses for Z\ = Z 2 vs. T. (a,b) n<i = ni/2; (c,d) n 2 = n\. The 
arrows indicate the mass average (m) = (mmj + n2Z 2 )/(nx + 712). The error bars represent 5% 
(10%) uncertainty in the measurement of the MD sound speed for M = 5 (20). 



C. Mode Dispersion 

Now we turn to the description of the full mode structure in the liquid state. By solving 
the characteristic equations for the matrices (f!2l) and ( Tl3l) one obtains the full u(k) dispersion 
for the 4 liquid modes. The results of this calculation are displayed for n 2 /nx = 1 and 1/2 
density ratios and for the already chosen Z = Z 2 /Zx = 0.7 and 1.4 (2.0), M = rriz/mx = 0.2 
and 5.0 parameter values. The Z 2 /Zx values have been chosen in the vicinity of the stability 
boundary for the (staggered rectangular and honeycomb) binary lattices. 

Our theoretical analysis of the mode structure was accompanied by detailed Molecular 
Dynamics studies of the dynamical fluctuation spectra of the system, as described below. 

In the Molecular Dynamics simulations we trace the trajectories of individual particles 

14 



as obtained from the integration of the their equations of motion: 



dv N 

where <pij is the interaction potential energy (oc ZiZj) between the particles % and j, and r/ij 
are the masses of the particles. We use periodic boundary conditions, the edge lengths of the 
computational box (L x and L y ) and the total number of particles are chosen to accommodate 
a perfect lattice for the selected density ratios (and expected associated lattice structures). 
In the case of n^jnx = 1 we use N\ = N2 = 2040 particles, while in the case of n^jny = 1/2 
we use Ni = 2720 and N2 = 1360 particles. 

In the simulations of liquid-phase systems normally random initial particle configurations 
are set up. In all cases ample time is given to the systems to reach thermodynamic equilib- 
rium before measurements on the systems start. During this equilibration phase, rescaling 
of the particle velocities is applied to reach the desired system temperature; this procedure 
is, however, stopped before data collection. 

The central quantities to be calculated in the simulations are the fluctuation spectra of 
the densities and currents. Static pair distribution functions gAB{ r ) = 1 + ^abO") are also 
obtained and used as input for the QLCA calculations. Information about the (thermally 
excited) collective modes is obtained from the Fourier analysis of the correlation spectra of 
the density fluctuations of the different species (A, B— 1,2): 

p A (k,t) = ^2exp[ikxj(t)], (31) 

yielding the dynamical structure functions as {59) ]: 

S AB {kM = —±^ ^^p A (k, U )^ B (k,^ (32) 

where At is the length of data recording period and p(k,u) = J-[p(k,t)\ is the Fourier 
transform of ( I3TI) . The (A,B) combinations label spectra related to component 1, Sn, to 
component 2, S 2 2, as well as to the cross term Si 2 . 

Similarly, the spectra of the longitudinal and transverse current fluctuations, L(k,u) and 
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T(k, oj) are obtained from Fourier analysis of the microscopic quantities, respectively, 

\ A (h,t) = ^2v jx (t)exp[ikxj(t)], 



3=1 



r A (k, t) = y~] v jy (t) exp [ikxj(t)] , 

3=1 



(33) 



where Xj and Vj are the position and velocity of the j-th particle. Here we assume that k 
is directed along the x axis. These calculations allow the determination of the spectra for a 
series of wave numbers, which are multiples of k min>x ^ = 2n/L x ^, where L x ^ is the edge 
length of the simulation box in the x (or y) direction. 

The identification of the collective modes is based on the observation of the extrema of 
Ln and L 2 2- When the peak positions do not completely coincide, (this may happen for 
various reasons, which will be discussed elsewhere), it is the position of the stronger peak 
that is accepted. 




FIG. 4. Liquid state: Examples for the g(r) pair distribution functions at T = 120 (a,c,e) ri2 = n\/2\ 
(b,d,f) ri2 = n\. Note that for Z\ = Z2 we obtain g\\ = #22 = 312, irrespective of the density ratios. 
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Distribution functions gABij) = 1 + }iab{t) that have been inputted in the QLCA calcu- 
lations are given in FIG H] for the previously chosen n^jnx and ZijZ\ values. We have also 
added the Z%jZx = 1 distribution functions, in order to show that in this case the three, 
hn, hi2, and h 2 2, correlation functions are identical, independently of the density ratios (the 
mass ratios obviously do not affect the correlation functions). 

Some illustrative current fluctuation spectra are given in Fig. 
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FIG. 5. Liquid state: Examples for the longitudinal current fluctuation spectra, (a) left column: 
ti2 = ni/2, Z2 = 2Zf, (b) right column n 2 = n\, Z2 = \AZ\. 



The MD simulated mode structures, together with the QLCA calculated dispersion curves 
are given in Fig. El Although the MD spectra are sometimes quite noisy as the collective 
modes have rather broad peaks in the spectra, the agreement between the simulated and 
calculated dispersions, in general, is good. A new feature shown by the simulation but not 
predicted by the QLCA formalism is the merging of a portion of the longitudinal acoustic 
and longitudinal optic modes at low m 2 /mi values into a new acoustic mode. 
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FIG. 6. Liquid state: Current fluctuation spectra from MD simulation (color map) compared with 
QLCA calculated dispersion (black lines) for T = 120. (a) ni = n\/2\ (b) ni = n\. 



ILL BINARY LATTICE 

Depending on the Z and n values of the two components, a variety of ordered and 
disordered phases should exist in a 2D binary crystal. In combination with the different 
melting temperature associated with the different phases, a rather complex phase diagram 
can emerge. The stability of the different structures can be analyzed through a thermody- 
l, 13] (minimizing the free energy) or through a dynamical normal mode 



namic approach 



IS 
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analysis. In this paper we restrict ourselves to the study of the T = (T — > oo) lattice 
structures only, which is amenable to the latter approach. 

The lattice calculation is based on the evaluation of the lattice sum for the dynamical 
matrix 



COT 



ZaZb 



j^OVij) (expHk • r hAB - Sab) 



S 



C^A j 



{ r j,AC> 

(34) 

over all the particle pairs with designated A, B and A, C indices, which now run over all the 
bases in the primitive cell (the number of which may be equal to or greater than the number 
of species, i.e. 2). The evaluation was done for ca. 10 5 particles. 

In the following we will consider two different lattice structures with the previously studied 
density ratios n 2 /n 1 = 1 and n%/n\ = 1/2. These two cases provide a reasonable guidance 
as to what lattice mode spectrum to expect in more general situations. In both cases we 
choose the equilibrium hexagonal lattice as the skeleton Bravais lattice. The descendent 
crystal structures should be stable in the vicinity of Z% = Z 2 . With Z\ — 1, Z 2 is restricted 
to Z m < Z 2 < Z M . The values of Z m and Z M have been determined by finding the onset 



of unstable normal modes 



60] and are given for both cases in Table [H Then the resulting 



lattice structures are the following; 

TABLE I. Stability regions. 



n 2 /m 


Zm 


Zm 


1 

1/2 


0.646 ± 0.001 
0.51 ±0.01 


1.548 ±0.002 
2.88 ±0.01 



1. In the equal-density case with rii = n 2 we obtain a staggered rectangular (SR) lattice 
with the aspect ratio : 1. 

2. In the half density case with n 2 = %/2 we obtain a honeycomb (HC) lattice for species 
1, while the particles of species 2 occupy the center sites of the honeycomb and form 
a hexagonal lattice; the lattice constants of the two lattices are in the ratio : 1, see 
figure [3 
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FIG. 7. Principal lattice structures: staggered rectangular (a), and honeycomb (b). Primitive cells 
are shown with dashed lines. In (b) the positions 1' and 1" are distinguished. 



The SR structure is built up from 2 bases in the primitive cell, while the HC structure 
has 3 bases in the primitive cell. According to 48|, |49[ other possible structures may exist 
outside the stability domains of Table [H such as various rhombic structures, the asymmetric 
hexagon (also with 3 bases in the primitive cell) and various pentagonal structures (with 3-5 
bases in the primitive cell). 



The number of modes r in general is r = d x b, where d is the dimensionality and b is 
the number of bases in the primitive cell. In general, the polarizations of the modes can be 
characterized only in the combined r-dimensional species-configuration space. In specific 
situations, however, (i) the r-dimensional space factorizes into the 6-dimensional species- 
and d- dimensional configuration sub-spaces; moreover, (ii) longitudinal and transverse po- 
larizations (with respect to k) may become the eigen-polarizations in the latter. This occurs 
when k is along one of the principal axes of the crystal. Thus the L + , etc. designations 
remain still meaningful and, by continuity, can be used for the labeling of the modes, with 
the proviso that since in general, more than one pair of optic modes may exist, a further 
index, say (3 = 1,11 may be needed for the full labeling. The HC mode structure consists of 
6 modes altogether, out of which 3 are "longitudinal" and 3 are "transverse" modes. Due 
to the rotational symmetry of the reciprocal lattice for k — > the L and T optic modes 
are degenerate at k = 0. In this limit, one can identify a pair of acoustic and two pairs of 
degenerate optic (gapped) modes. The SR mode structure consists of 4 (2 longitudinal and 
2 transverse) modes. In the absence of rotational symmetry the L and T modes are not 
degenerate at k = 0, and the lo^_ and uJ^ gaps are separated. 
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A. Optic modes 



The simple geometric structure of the primitive cell allows one to obtain a transparent 
result for the u(k — Y 0) frequency gaps. The results are given below and portrayed in Fig. 
[TU1 For the HC at k = the elements of the C-matrix are 

CxV(O) 



C*,(0) = J_ 



, < 2 
u 1 



cm 



"2 



C^(0) 



1 z 

2V2m 
1 

2v/2 
Z 



2V2m 



^^(r i)21 ,)+* x (r i|21 „) 
i 

S^fo.i'i" 



(35) 



7M'2 



<7& 2 (0) 



By symmetry, all the lattice sums are equal; the rotational symmetry (L = T) can be 
further exploited to obtain 



y 



J.21 

in terms of which the roots of the cubic equation are 



(36) 



L.T 



T j = v/2(p 2 + g 2 )P 



(37) 



Note that is an "invariant mode", where the gap frequency is independent of m 2 ; 

in this mode the light particles oscillate around the inert heavy particle. 
For the SR a similar construction yields 



Q ' 



L,T 



1 



2^2 



.i,i2 



(38) 



in terms of which 



(39) 



-T 



Here the P, Q-s are lattice sums, characteristic of the lattice structure (SR or HC); they 
depend on k only. They can be contrasted with F 12 factor appearing in the gap frequency 
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expression in the liquid ffTB"]) , that depends on Z 2 jZ\ as well (see Fig. [9]), but for Y — > rf reeze 
its value in the n 2 = \j2n\ and n 2 = n\ cases reasonably well approaches the corresponding 
4P (for the HC) and 2(Q L + Q T )/2 (for the SR) values respectively. 

While the gap frequencies are angle independent, the polarizations associated with them 
are not: Fig. [8] shows that T and L polarizations switch place as the propagation angle 
varies from to 90 degrees. 
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FIG. 8. SR lattice: polarizations of the gap frequencies versus propagation angle for Z 2 = Z\. 



Figure |9] shows the Z and m dependences of the respective gap frequencies in the SR and 
HC crystal lattices and the corresponding gap frequency in the liquid. Commenting on the 
HC case first, we note that the liquid has only one frequency gap and therefore there is no 
equivalent of the invariant mode in the liquid. Turning to the SR lattice, one observes the 
separation of the longitudinal cj^. and transverse gaps. The Co>gap frequency in the liquid 
largely follows the angular average of U\ and u 2 , but less closely than it does in the case of 
the YOCP [61 1 . 



Figure [10] shows the dependence of the P, Q lattice sums on the screening parameter; 
the smooth extrapolation to the k = value provides the input for the calculation of the 
noteworthy Coulomb gap frequencies via Eqs. [37]and|39l In parenthesis we note that a little 
reflection shows that the P(k = 0) value bears a close relationship to the M = r~ 3 dipole 



sum over a hexagonal lattice whose value is well-known 62]: M/2 = 0.7985/6 3 in terms of 



the Wigner-Seitz radius b. Then P(k = 0) = 2- 13 / 4 (3 3 / 2 - l)M/2. 
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FIG. 9. Mass and charge ratio dependence of the gap frequencies. The QLCA and MD result are 
also shown, (a) HC lattice; note the portrayal of the "invariant mode" in the right (Z = const.) 
panels, (b) SR lattice. 



B. Acoustic modes and sound speed 



In contrast to the optic modes, whose dispersion is highly structure dependent and is, 
in general, quite different from the corresponding mode dispersion in the liquid, the k — > 
behavior of the acoustic phonons in the lattice is largely similar to that of their liquid 
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FIG. 10. The dependence of the lattice sums on the Yukawa screening parameter (k). (a) HC 
lattice; (b) SR lattice. 



counterpart. More precisely, the sound speeds, as calculated by the QLCA and verified by 
simulations, go over quite smoothly to the lattice sound speeds as T crosses the freezing 
boundary. This is visible in Fig. [2j The only difference of some significance arises in the 
case of the SR lattice, due to the fact that its reciprocal lattice space is, in contrast to the HC 
structure, anisotropic even in the k — > limit. The most important observation, however, 
is that the notion of the VAA as a dominant feature for the low frequency excitations both 
in the liquid and in the solid state, is of universal validity. 



C. Mode dispersion 

The full calculated lattice phonon dispersion diagrams both for the HC and SR lattices 
are portrayed in Fig. [TTJ In order to be able to compare the MD results with lattice 
summation data, simulations were carried out at very low temperatures, at Ti = 10 4 . In 
these runs the particles are initially arranged in a perfect lattice and their thermal motion 
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does not disrupt the lattice in the course of the simulations. In Fig. PT21 we display the MD 
simulation results for these finite temperature lattices: the MD simulations and the results 
of the lattice calculations are in full agreement. 

dog 90 deg 




FIG. 11. Calculated mode dispersions for different propagation angles, (a) HC lattice, note that 
the invariant mode frequency at k = (pointed at by the arrows) remains invariant for any 
M = milm\ and any angle; (b) SR lattice. 



The polarizations of the modes in the combined (in general not factorizable) cartesian and 
species space can be assessed from Figures [131 [TH and CLS, [T6], where the components of the 
eigenvectors for the HC (SR) lattice modes along the 6 (4) eigendirections of the dynamical 
matrix for a given k are shown. The lengths of the La and T4 labeled bars (components 
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of the eigenvectors) are proportional to the longitudinal and transverse displacements of 
particles at position A. Samples are given for propagation angles along and off the principal 
axes. Note that in the latter case no overall polarization direction can be assigned to the 
displacement of the particles belonging to different species. 

Finally we address the question of how the collective mode dispersion depends on k, 
the screening parameter of the Yukawa potential and, in particular, how the transition to 
the k = Coulomb limit occurs. Figure [T7] shows that there is a smooth evolution of the 
mode dispersions towards the Coulomb limit and towards the changeover of the longitudinal 
acoustic mode into the characteristic quasiacoustic \/k Coulombic behavior. It will be shown 
in another publication, that this behavior is in sharp contrast to what happens in the 3D 
case [ss]]). 



IV. COMPARISONS AND CONCLUSIONS 



The results of earlier analyses 1CH15| of the collective mode structure of the YOCP have 
established the close affinity of the mode structures in the strongly coupled liquid and in 
the crystalline lattice states. More precisely, what has been found is that the QLCA model, 
which essentially portrays the strongly coupled liquid as a superposition of randomly oriented 
microcrystals and determines the eigenmodes as those of the averaged crystal, provides an 
adequate description of wave propagation in the liquid. Whether such a simple picture would 
prevail in the binary liquid where the non-random distribution of the particles belonging to 
the two species is an issue as well, is not a priori obvious. The study presented in the previous 
Sections shows, however, that this is the case. In the following, we discuss the relationship 
between the phonon dispersion in the binary crystal, as calculated by lattice summation and 
corroborated by MD simulations, and collective excitations in the binary liquid, as provided 
by the QLCA description and the MD simulations. Judged by comparison with the results 
of the MD simulations, the QLCA results are quite reliable, with two exceptions, that we 
will discuss below. In comparing mode structures in the liquid and in the solid, the effect 
of the different density ratios has to be kept in mind: while in the former the difference 
between the HI2 = n\ and ri2 = na/2 cases does not make a major difference, in the latter 
the two different crystal structures (SR and HC) substantially affect the mode structure. 

Focusing first on the low-/c acoustic excitations, we see that there is an almost perfect 
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FIG. 12. Current fluctuation spectra from MD simulation at Ti = 10,000 (color map) compared 
with dispersion from lattice calculations (black lines), (a) HC lattice; (b) SR lattice. 



agreement between the liquid QLCA and MD sound speeds, on the one hand, and the 
corresponding values in the liquid and in the two crystal structures studied, on the other 
hand. The only difference of note, as we have already pointed out, arises in the case of 
the SR lattice, due to its anisotropy, which results in a narrow band of sound speeds; in 
the liquid, as represented by the QLCA, it is replaced by an angular average. The most 
important result that emerges from all this is the fact that the low frequency excitations 
are governed by oscillation frequency u of the Virtual Average Atom (see Eq. 126]) which 
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FIG. 13. HC lattice: Mode polarizations for deg propagation. Z = 0.7, M = 5, a = 0° (A;||x), 
ka = 0.2 in the panels, particle "2" is the heavy one (see Fig. [7]). 



is created by the average charges and masses of all the components. This effect, as it was 
discussed in some detail elsewhere jsij , has its most dramatic manifestation for the effective 
mass of the nominal plasma frequency of the binary (with respective masses mi and r«2), 
which in the weakly coupled case is formed, in general, through the "parallel connection" 
of the two masses (l/m c g = \jm\ + l/m 2 ), but which in the strongly coupled case becomes 
the "series connection" of the two masses (m e ff = vn,\ + m 2 ). While the VAA has been a 



useful heuristic concept for liquid alloys 64j and for disordered systems 65H68], and also in 
connection with self-diffusion 53|, here we have been able to give a rigorous demonstration 
through the QLCA of the emergence of this phenomenon. The MD simulation has shown (see 
Figs. [T] [3]) that in the T — > r freeze limit the VAA concept becomes "exact", in the sense that 
after the subtraction of the explicitly identifiable pair correlation /i 12 dependent correlational 
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contribution it determines the sound speed. With decreasing T the m e g decreases, seemingly 
marching towards weak coupling limit, but within the boundaries of our MD simulation 
which covers only the T > 5 domain, the behavior is still essentially strongly coupled, in 
that the decrease of m e s from its high T value is quite slow. However, this decrease of m e s 
in the moderately coupled domain is not reflected by the QLCA model: there m e g preserves 
its high T value (Eq. 129]) for any T. This is the first inadequacy of the QLCA and it is the 
consequence of the fact that the appearance of the VAA structure is formally correlation 
independent. Correlational effects appear only indirectly, through the model from which it is 
derived and which adopts quasilocalization as its basis. That the quasilocalization can lead 
to such a qualitative effect is a novel feature of the approximation, which manifests itself 
only in binary systems. In contrast, in the single component system, the weakly coupled and 
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FIG. 15. SR lattice: Mode polarizations for deg propagation. Z = 0.7, M = 5, a = 0° 

ka = 0.2 in the panels, particle "2" belongs to species "2". Note the L/T and 1/2 polarization 

mixings. 

strongly coupled states differ through their explicit correlation function dependence only. 

A hallmark of the binary system is the emergence of - one or more - optic modes with 
a k = gap frequency. In the liquid state there is only one gap frequency, corresponding 
to the two - longitudinal and transverse - modes that become degenerate at k — 0, due 
to the isotropy of the liquid. In the crystal lattice this degeneracy may or may not be 
lifted, depending on the local environment: it is in the SR crystal, but it survives in the 
HC crystal. In addition, in the crystal lattice the number of optic modes increases with 
the number of particles in the unit cell, which increases in oder to accommodate n 2 /n 1 
unequal 1 density ratios: hence an additional degenerate gap frequency in the HC crystal. 
This latter is the "invariant mode" whose gap frequency is independent of the mass of the 
lower density component, which remains inert in this mode. The mode does not have an 
equivalent in the liquid. The other, "normal" mode does re-appear in the liquid, with the gap 
frequency in the vicinity of the crystal equivalent (for the HC) or between the longitudinal 
and transverse gaps (in the SR). It should be emphasized though that the approximation of 
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FIG. 16. SR lattice: Mode polarizations for 10 deg propagation. Z = 0.7, M = 5, a = 10°, 
ka = 0.5 in the panels, particle "2" belongs to species "2". 

the liquid dispersion by angle averaging the lattice phonons is not equivalent to the QLCA. 
This difference was already demonstrated for the YOCP: here it is much more pronounced. 

The gap frequencies are not related to the VAA. In the liquid they can be expressed in 
terms of the nominal Einstein frequencies VLab (Eq. 115]) and thus they follow the "parallel 
connection" rule. 

In the liquid state one can identify two upper (Clj) and lower (flu) Einstein frequencies. 
(Eq. [T9|) For high k values the two "acoustic" (longitudinal and transverse) modes of the liq- 
uid merge into Clu, while the two optic modes merge into Clj. These latter cannot be directly 
identified in the crystal lattice, but they appear in the expressions for its gap frequencies, 
showing a good agreement with the QLCA calculated liquid Cli and tin quantities. 

According to the MD simulation result (Fig. [6]) the slopes in the vicinity of k = 
of the longitudinal acoustic and longitudinal optic modes match and the two modes fuse 
into a single acoustic mode. There is no indication of this phenomenon within the QLCA 
formalism. 

As to the dependence on the screening constant k, we see (Fig. [17]) that the qualitative 
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FIG. 17. The dependence of the HC lattice dispersions on the Yukawa screening parameter (k) at 
m,2 = 5mi and Z2 = 0.7 Z\. 

features of the dispersion remain unaffected over a wide range of k values, down to and 
including the k — Coulomb limit. 

A number of problems relating to the collective dynamics of the system have been iden- 
tified, but have not been studied in this paper: the damping of the modes, the detailed 
structures and the link between the various fluctuation spectra, the nature of the underlying 
order in the liquid phase, lattice stability and structures, etc. These problems will have to 
be investigated in future works. 
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